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The fundamental difficulties associated -with ab initio 
solutions of the quantum mechanical many-body prob- 
lem stem directly from the large dimensionality of the 
-wave function. Rather than dealing directly -with the 
many-body -wave function, large scale electronic struc- 
ture calculations determine the ground-state energy by 
minimizing an energy functional over the much more 
manageable space of single-particle orbitals. Hohenberg- 
Kohn-Sham theory pp] places this approach on a firm 
theoretical footing by establishing that exact electronic 
ground-state energies can be determined in this manner 
by minimizing a universal functional of such orbitals. 

Although the requisite universal energy functional is 
not kno-wn exactly, the local density approximation |^] 
and various improvements thereon |^-|^ are sufficient to 
resolve bond energies to a little better than one-tenth of 
an electron Volt. With this accuracy, these functionals 
make possible highly predictive first principles studies in 
diverse areas such as the study of surfaces, point defects, 
plastic deformation and chemical reactions. (For a revie-w 
see 1^.) Although sufficient for many studies, this error is 
too large to allo-w accurate prediction of the rates of mi- 
croscopic processes at room temperature, thus limiting 
the ultimate predictive po-wer of the ab initio density- 
functional approach. Whereas further improvement of 
energy functionals for single-particle theories is an active 
area of research, much less -work has been done to con- 
struct energy functionals of the one-body density matrix 
1^,^ . Such functionals arc promising because the density 
matrix contains more explicit information than do the 
Kohn-Sham orbitals, and an accurate energy functional 
thus is likely to be simpler in form. The exact kinetic 
energy functional, for instance, is kno-wn for the density 
matrix, but not the Kohn-Sham orbitals. 

To date, the only extensively used density-matrix en- 
ergy functional is the Hartree-Fock approximation. One 
may vie-w Hartree-Fock theory as approximating the f-wo- 
body matrix as a sum of Hartree and exchange terms, 
each of which is a tensor product of the one-body den- 
sity matrix. Hartree-Fock then uses the known, exact 
energy functional of the two-body density matrix to pro- 



duce an energy functional of the one-body density ma- 
trix. In this work, we build upon this perspective and 
consider tensor product approximations to the two-body 
density matrix which not only serve as generators of en- 
ergy functionals of the one-body density matrix but also 
provide estimates of the two-body matrix in terms of the 
one-body matrix and thereby shed light on the nature of 
correlation in quantum many-body systems. 

Notation — For a system of N electrons, the ex- 
act energy functional of the two-body density matrix 
7(xiX2,xix^) = (vI/|V't(xi)V't(x!>)^(x2)V;(xi)|vI/)/2, is 
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Here and throughout, we work in atomic units, and the 
Xi are compound coordinates representing both position 
Ti and spin so that integration over x^ represents in- 
tegration over space and summation over spin channels. 
The external potential and one-body density matrix are 
C/(r) and n(xi,x'i) = (^'|i/)'l'(x'i)Vi(xi)|5'), respectively. 
Finally, the one-body density matrix comes directly from 
7 through the sum-rule 
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n(xi,x'J = ydx2 7(xiX2,x'iX2) 
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Computational considerations — Although is ex- 
act, the two-body density matrix 7 is a function of four 
variables, and direct computation with such functions is 
infeasible. Viable techniques, however, exist for deal- 
ing directly with two-variable functions such as the one- 
body density matrix [^-^|. Accordingly, we expand 7 
in terms of two-variable functions, 
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which, by separation of variables, is always possible. 
Here, gi and hi are one-body operators (functions of two 
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variables), we view the two-body density matrix as a 
two-body operator (a function of four variables) , and the 
tensor product denotes one of the three possible choices 
for separating four variables, 

g<»h^ 5r(xi,Xi)/l(x2,X2) 

gf^h^ 5(xi,X2)/i(x^,x'i) 
g<»h = .g(xi,x^)/i(x2,x'i). 

Quantum statistical considerations — We wish to max- 
imize the physical content of a truncation of the ex- 
pansion to a limited number of terms. Physically, 
n(xi,x'i) is the quantum amplitude for the insertion of 
a hole at position xi and its instantaneous removal from 
x']^, whereas 7(xiX2, x'j^Xj) is the quantum amplitude for 
the insertion of a pair of holes at Xi , X2 and their removal 
from x'j^jXj. Approximating the latter events as inde- 
pendent gives 7 « jH = {h®h)/2, the familiar Hartree 
approximation, where the factor of two maintains the 
normalization implicit in (^. We may then refine this 
mean-field behavior and define 7 = + 7xc, where 
represents exchange and correlation effects which we then 
expand, without loss of generality, according to 

The most significant drawback of the Hartree approx- 
imation is that 7// is not properly antisymmetric with 
respect to particle exchange (x'j^ ^ Xj), so that, poten- 
tially, many terms will be required in to restore this 
symmetry. Alternately, we may explicitly ensure anti- 
symmetry and take 7 = ^hf+1c = {n®n—h®Ji) 
In this form, ^hf is precisely the familiar Hartree-Fock 
approximation. Again, we may expand the unknown 7c 
as in (||). 

Finally, we note that a fundamental difficulty exists 
when using a finite number of terms of Type II. Because 
of the instantaneous nature of the quantum event which 
7 represents, we expect for normal systems that 7 — > 
as the insertion x']^X2 and removal X1X2 locations are 
placed at ever further distances from one another. In 
an extended system with the expansion for 7 containing 
a finite sum of terms of Type II, taking this limit while 
keeping the insertion points near one another and the 
removal points near one another will violate this asymp- 
totic condition unless each Type II term is identically 
zero. Accordingly, we focus below on approximating 
and 7c with terms of Type I and III. 

Symmetry — Physical symmetries of 7 significantly re- 
strict the allowable forms for the tensor products appear- 
ing in expansions of the form (j^) for 7xc and 7c. We first 
consider Hermiticity 7(xiX2,X2X2) — 7*(x'^X2, X1X2) 
and particle permutation symmetry 7(xiX2,X2X2) = 
7(x2Xi, Xjx']^) for each of the three types of tensor prod- 
uct and defer discussion of Fermionic antisymmetry until 
after the sum-rule immediately below. Because both 7// 
and ^HF respect these symmetries, so must tensor prod- 
ucts representing and 7c. These symmetries imply 
relations with arbitrary constants of proportionality be- 
tween the one-body operators of such products which 



may be absorbed into the one-body operators. Table || 
summarizes the final result of such considerations for all 
three tensor products. We note that in all cases symme- 
try restricts the tensor product to involve only a single 
one-body operator which remains to be determined. 

Sum-rule — Remarkably, when a single term satisfy- 
ing Hermiticity and particle permutation is used to rep- 
resent 7a;c or 7c, we can "invert" the sum-rule and 
determine the two-body density matrix 7 in terms of the 
one-body matrix n. Inserting the resulting form for 7 
into (|l|) then generates an energy functional of the one- 
body matrix which implicitly satisfies the sum-rule, an 
important property for obtaining good ground state en- 
ergies. 

When terms of Type I represent ^xc or 7c, the sum rule 
becomes ^(Tr^) = 6, where = n or n(l — n), respec- 
tively. The solution of this equation \s g — ±6/vTr6. 
In either case, the extensivity of Tr 6 makes the resulting 
solution irrelevant in the thermodynamic limit N ^ 00. 
For products of Type III, the sum-rule combined with 
appropriate symmetries gives if' = 6, so that g ~ \fd 
with o defined as above. 

Fermionic antisymmetry — Fermionic antisymmetry, 
7(xiX2, x']^X2) = — 7(xiX2, Xjx']^), is a stronger condition 
than the particle permutation symmetry condition con- 
sidered above and imposes more complex constraints. 
Two-body density matrices 7 which satisfy this condi- 
tion may always be written with Type I and HI products 
appearing in corresponding pairs of the form g®g — g®g 
and with Type II products (were we to consider such) 
whose one-body terms are separately antisymmetric. 

The two viable functionals remaining after the above 
considerations involve a Type III representation for either 
Ixc or 7c. Unfortunately these forms do not consist of 
symmetric pairs of Type I and Type HI products. To sat- 
isfy antisymmetry, one could take 7a;c = —{h®n)/2, so 
that 7 = {n®n — n®n) /2, which is simply the Hartree- 
Fock approximation. To go beyond this, antisymme- 
try requires representing 7c as a pair of terms so that 
7 = {n®h — h®n + g®g — g®^g)/2. 

The sum-rule for this latter extension on Hartree-Fock 
is if' — g(Tr g) = h(l — h), which has no solution for g in 
extended systems unless Tr g vanishes in the thermody- 
namic limit. Under this condition, we have the solution 
g — ±-^ri(l — n), where the signs of the eigenvalues in 
the square root must be chosen to ensure Tr 5 ~ 0. For 
the paramagnetic phase we consider below, the natural 
choice is to take the upper sign for the spin-up block 
of the density matrix and the lower sign for the spin- 
down block. The structure of the energy functional (|^), 
however, is such that for this choice the resulting energy 
functional of n is equal to that generated by representing 
7c as a single Type HI product. 

Representation of the functionals — The above dis- 
cussion leaves the energy of two representations of the 
two-body matrix to explore: corrected Hartree theory, 
7 = IH — {Vfi®Vn)/2, and corrected Hartree-Fock the- 
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ory, 7 = ^HF — — n)^^^h{\ — h))/2. The cor- 

responding energy functionals may be represented either 
directly in terms of the one-body density matrix for use 
with direct density-matrix methods ]lO[-p^, or in terms 
of the spectral ( "natural orbital" ) representation of the 
density matrix. 

The energy functional contains one-body terms 
which may be evaluated directly in terms of the one- 
body density matrix. The remaining two-body term 
may be computed in terms of the two-point density 
c(ri,r2) = X]siS2 ■^(■''■i-''-2; ^1^2)- When the above func- 
tionals represent paramagnetic states, this density has 
the following form in terms of the eigenvectors (natural 
orbitals) 0i and corresponding eigenvalues (occupancies) 
/,; of the one-body density matrix fi, 



c(ri,r2) 



4E/^/j'^^('-i)<^*(^i)'^^-(^2)0*(r2) (4) 



where for corrected Hartree and Hartree-Fock the- 
ories, we take = \/ hf] and g^j = /j/^ -|- 
- - /j), respectively. 

The first of the above forms corresponds to the natural 
orbital functional proposed by Goedecker and Umrigar 
1^ as an ansatz among the many possible choices for the 
Qij which satisfy the sum-rule. Here, we see that this 
form is the unique Type III correction to Hartree theory 
which satisfies basic symmetries and the sum-rule. This 
theory has been studied analytically in the case of the 
homo gen eous electron gas in the low density regime (r^ > 
5.77) iQ. Here, we present results for both this theory 
and our new corrected Hartree-Fock theory over a range 
of densities which includes both this low density regime 
and the regime more physically relevant in the valence 
region of solids, rs < 5.77. 

Homogeneous electron gas — Translational invariance 
in the homogeneous gas ensures that the natural orbitals 
are plane waves. The remaining freedom in h lies in the 
occupancy eigenvalues /k. For the paramagnetic phase, 
the energy then takes the following form. 



V'^dkdk' 



-ffkk' 



(5) 



where V is the volume of the system and the factor of 
two in the kinetic energy arises from the sum over spin. 

The minimum of this energy functional occurs when 
S{E — fiN)/6f-k = 0. For both corrected Hartree and 
corrected Hartree-Fock, this gives to leading order as 
|fc| 00, 

d{E-fiN)/6h^k^-C/^, 

where C is a constant. Hence, for large fc, the momen- 
tum density scales as /k ~ |k|"^, in agreement with the 



random phase approximation and other many-body cal- 
culations This analysis holds for the present two 
functionals at all particle number densities and, in fact, 
for all theories of the form of (^) in which gkk' ^ \//k./k' 
as fc — > 00 and /k — > 0. 

Numerical Results — To determine the total energy 
of the homogeneous electron gas within the corrected 
Hartree and Hartree-Fock approximations, we exploit 
spherical symmetry and reduce /k to a single variable 
function /(fc), which we represent on a radial mesh of 
variable spacing. We then minimize (|^) with conjugate- 
gradients techniques subject to both the constraint that 
f{k) corresponds to a total of N electrons and the Fermi 
constraint that < f{k) < 1. The latter constraint we 
implement by defining f{k) = (e~^ 4- 1)~^, where x is 
then free to range over the real line. 

To evaluate the energy functional in terms of the fill- 
ings, we integrate numerically the kinetic and exchange- 
correlation terms in (|^). A key step in producing ac- 
curate results is to avoid numerical integration across 
the familiar logarithmic singularity at k — k'. For cor- 
rected Hartree theory, we avoided this by integrating 
the exchange-correlation terms once by parts analytically 
and then evaluating the result numerically. For the cor- 
rected Hartree-Fock case, we found it more efficient to 
approximate the Coulomb potential by a Yukawa poten- 
tial, e~'"'/r, and numerically extrapolate the results to 
K = 0. For these latter calculations, we found it more effi- 
cacious to work in the grand canonical ensemble (holding 
the chemical potential fixed while minimizing E — fiN) 
than to impose the particle-number constraint explicitly. 

The above procedures introduce four sources of numer- 
ical error: the finite spacing of the radial mesh, the finite 
extent of the mesh in fc-space, the error in the extrapola- 
tion in the case of using the Yukawa potential, and the fi- 
nite number of iterations involved in our searches for the 
minima. Through analysis and numerical experiments, 
we determine these errors to be less than 10"'^, 10"*^, 
10"'^ and 10~^ Hartree per particle, respectively, over the 
range of densities which we explored (0.01 < rg < 30). 

We found the minimization of the above functionals 
problematic for two reasons. First and foremost, the dis- 
cretized representation of the one-dimensional function 
f{k) on the radial mesh artificially stabilizes any discon- 
tinuities in f{k) which may arise during the minimiza- 
tion process. This may be managed by taking care to 
ensure that the radial mesh provides sufficient resolution 
to describe any large gradients which arise in f{k). An 
additional factor making the calculations difficult is that 
both the corrected Hartree and corrected Hartree-Fock 
functionals suffer from a very large condition number (in 
excess of 10^). 
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FIG. 1. Correlation energy per particle of the uniform 
electron gas. Hartree-Fock (horizontal dashed line); exact 
many-body results 1^,^ (dashed curve), corrected Hartree 
theory (open circles), corrected Hartree-Fock theory (dia- 
monds), analytic results for corrected Hartree theory in 
the low density limit (solid curve). 

Figure |^ compares our results for the correlation en- 
ergy per particle of the electron gas with established 
many-body predictions. We find that at very low den- 
sities, both new theories begin to approach one another 
as — + oo. For even moderately low densities (r^ > 6), 
both theories perform significantly worse than Hartree- 
Fock, overestimating the correlation energy by factors 
from two to four and higher. For the densities experi- 
enced in the valence regions of solids (6 > Ts > 1), cor- 
rected Hartree-Fock theory (diamonds) crosses over from 
over-correlation to under-correlation, whereas corrected 
Hartree theory (circles) over-correlates by roughly the 
same amount by which traditional Hartree-Fock theory 
under-correlates. For higher densities more representa- 
tive of atomic cores (1 > r^), the corrected Hartree- 
Fock theory predicts near zero correlation. At these 
same densities, corrected Hartree theory remains over- 
correlated and gives a somewhat better estimate of the 
correlation energy than does traditional Hartree-Fock. 
The asymptotic form of the corrected Hartree results in 
the high density limit appears to be a straight line in our 
plot, which corresponds to the form of the leading order 
terms of the Gell-Mann-Brueckner expansion . How- 
ever, the prefactors which we find by a least squares fit, 
E « 0.0570 log(rs) — 0.0714, are only correct to within 
a little better than a factor of two: the well-known ana- 
lytic values for these coefficients are 0.0311 and —0.0480, 
respectively. This suggests that the favorable reports for 
atoms for the natural orbital functional |3j may be in 
part a consequence of the high densities experienced in 
atomic cores, the limited particle number in the low den- 
sity regions, and the self-interaction corrected nature of 
those calculations. 
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FIG. 2. Predictions of momentum distribution of the uni- 
form electron gas for = 7 and Vs — 1: corrected Hartree 
theory (dashed curves), corrected Hartree-Fock theory (dot- 
ted curves), analytic results for corrected Hartree theory 
in the low density limit (solid curve). Inset shows magnified 
view. 

Figure ^ shows momentum distributions for = 7 
and rs ~ 1, which are representative of what we find at 
low and high densities, respectively. We note that the 
numerical results of both theories exhibit the correct 
scaling at large momenta and that we have very good 
agreement with the analytic results for the natural orbital 
functional (applicable for Tg > 5.77). 

We observe that at low densities the momentum distri- 
butions of both corrected Hartree theory and corrected 
Hartree-Fock theory resemble one another quite closely 
and, in fact, approach one another as rg approaches in- 
finity. This occurs because in this limit both theories 
predict very low occupations fi << 1 so that the cen- 
tral quantities in the two theories approach one another, 

The next panel of the figure illustrates that for densi- 
ties outside the range of the applicability of the analytic 
solution {rs < 5.77), this solution predicts / > 1 for suffi- 
ciently low momenta. We find that for such densities, the 
momentum distributions of corrected Hartree and cor- 
rected Hartree-Fock theory exhibit different behaviors at 
low momenta. Our numerical results indicate (to bet- 
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ter than 1 part in 10^) that the momentum distribution 
of corrected Hartree theory saturates to / = 1 near the 
Fermi level. (See inset.) Such behavior is in direct con- 
tradiction of analytic and numerical many-body results 
for the electron gas ||l^,|l8| and is contrary to the previ- 
ously reported experience with atoms for this functional 
1^. In contrast, the momentum distribution of corrected 
Hartree-Fock theory never saturates the Fermi occupa- 
tion constraint, but rather approaches a value at fc = 
which remains less than unity. We find that in the high 
density limit the corrected Hartree-Fock momentum dis- 
tribution in fact approaches the standard Hartree-Fock 
result, ultimately leading to zero correlation energy for 
this theory in the high density limit. 

In conclusion, we find that when the two-body den- 
sity matrix is expanded beyond traditional Hartree and 
Hartree-Fock theories in terms of tensor products of 
one-body operators, fundamental symmetries allow the 
exchange-correlation hole sum-rule to be inverted to give 
estimates for the two-body matrix, and thus correlations, 
in terms of the one-body matrix. By approximating the 
two-body density matrix in terms of the one-body den- 
sity matrix, this approach also leads to new energy func- 
tionals of the one-body density matrix. 

The addition of the first separable term beyond 
Hartree theory shows some promising features when ap- 
plied to the homogeneous electron gas, such as giving 
the correct scaling behavior of the momentum distribu- 
tion in the high momentum limit and the correct scaling 
behavior of the correlation energy in the high density 
limit. However, the electronic states predicted by this 
functional are significantly over-correlated, and attempts 
to improve upon it by the addition of further terms with- 
out additional constraints would only add more freedom 
in the minimization and lead to further over-correlation. 
Hence, before adding a second term, it is important to 
use antisymmetry to constrain the first such term to be 
the traditional Fock term. This gives the second func- 
tional which we considered, corrected Hartree-Fock the- 
ory. This second theory under-correlates the electron 
gas in the high density limit but still over-correlates at 
lower densities. This indicates that the addition of fur- 
ther terms will only hold the promise of producing a more 
accurate functional if additional, appropriate constraints 
are imposed. 
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Hermiticity 


Particle Permutation 


g®h 


g ^ g^h^h'' 


g = h 


g®h 


g = h^ 


g = ±f, h = ±h^ 


g®h 
III 


g = ±/V 


g ^h 



TABLE I. Implications of Hermiticity and particle permu- 
tation symmetries for tensor product expansions of 7. 



ACKNOWLEDGMENTS 



CSG would like to thank Sohrab Ismail-Beigi for very 
useful discussions. This work was supported primarily by 
the MRSEC Program of the National Science Foundation 
under award number DMR 94-00334. 



5 



